Growth exponents reflect evolutionary processes and treatment response in brain metastases

Tumor growth is the result of the interplay of complex biological processes in huge numbers of individual cells living in changing environments. Effective simple mathematical laws have been shown to describe tumor growth in vitro, or simple animal models with bounded-growth dynamics accurately. However, results for the growth of human cancers in patients are scarce. Our study mined a large dataset of 1133 brain metastases (BMs) with longitudinal imaging follow-up to find growth laws for untreated BMs and recurrent treated BMs. Untreated BMs showed high growth exponents, most likely related to the underlying evolutionary dynamics, with experimental tumors in mice resembling accurately the disease. Recurrent BMs growth exponents were smaller, most probably due to a reduction in tumor heterogeneity after treatment, which may limit the tumor evolutionary capabilities. In silico simulations using a stochastic discrete mesoscopic model with basic evolutionary dynamics led to results in line with the observed data.


INTRODUCTION
Macroscopic tumor growth is a complex process resulting from the interplay of different biological elements at the cellular and subcellular levels. These include the driving molecular alterations and their associated heterogeneity, angiogenesis, the immune system, the tumor microenvironment, and surrounding healthy structures, the effect of treatments on the different tumor phenotypes/genotypes, etc.
Mathematical growth laws have been shown to describe longitudinal tumor growth dynamics effectively in simple experimental models [1][2][3][4][5] . A great deal of data is available for those models, which do not have the biological complexity of human tumors. One would expect that describing cancer growth mathematically in humans would be far more difficult because of the different biological mechanisms that drive it over distinct tumor stages.
Assessing complex tumor growth dynamics over time, given the large genotypic and phenotypic variability developed during tumorigenesis, is very difficult with current techniques. Medical images are performed routinely in most cancer patients and provide rough global macroscopic information -the so-called imaging phenotype-which integrates the several processes occurring at the microscale, potentially providing information on the underlying tumor biology. However, longitudinal datasets for untreated malignant tumors are rare, and of limited quality, since therapeutic action is typically performed promptly. This is why studies of untreated tumor dynamics in humans have been mostly limited to the use of two time points per patient which allows only for the identification of average growth rates 6,7 .
A recent study using more time points has shown that the explosive growth observed in human cancers may be the result of evolutionary dynamics, which underpins increased tumor aggressiveness and leads to an acceleration in tumor growth as the disease progresses 8 . It is well known that solid tumors contain subpopulations of clonal cells with different genomic alterations, a phenomenon referred to as intra-tumor heterogeneity 9 . A growing body of evidence has highlighted the importance of understanding tumor heterogeneity and its evolution in cancer research 10,11 . According to this recent work 8 , the sequential emergence of new, more aggressive clones would produce a sustained increase in the average growth exponent of the entire tumor, as these clones would eventually prevail over the rest by competitive advantage, increasing their relative abundance and imposing their faster growth. This sustained generation of variation within a tumor has been shown to play a crucial role in the development of many aggressive clinical characteristics of cancer 12 . As such, measuring tumor evolution and heterogeneity can have significant clinical implications for cancer therapy 13 .
Our goal in this study combining clinical data and mathematical models, was to provide more evidences on the macroscopic growth dynamics of both treated and untreated tumors in human patients. We focused on growth dynamics of brain metastases (BMs) integrating mathematical and computational models with data analysis methods on a patient database with longitudinal magnetic resonance imaging (MRI) follow-up and more than a thousand lesions. BMs are the most common intracranial tumor and a major complication of many cancers, with 20-30% of cancer patients developing the condition in the course of their disease 14,15 . Here we focused on two scenarios. The first one was the growth of untreated BMs. Human beings as well as mice data were studied. The second scenario considered BMs under different treatment modalities. Therapy options for BMs are surgery, whole brain radiation therapy (WBRT), gamma knife radiosurgery, stereotactic radiosurgery (SRS), targeted and systemic therapies.
We wanted to investigate whether there was a measurable effect on the growth dynamics of the expected loss of tumor heterogeneity after treatment. To this end, we performed an analysis of the growth dynamics of untreated BMs, and BMs treated under different types of therapy, by fitting a general tumor growth law 16 . The results obtained were endorsed by in silico simulations with a discrete stochastic model of tumor growth, which helped to gain insight into the causes of the behavior observed in BMs.

Brain metastasis growth is explosive in untreated patients but not in treated ones
The Von-Bertalanffy model 16 , has been proven recently to describe the longitudinal dynamics of growing human tumors 8 . Eq. (1) describes the change in size (i.e., volume change dV dt ) experienced by an organism according to its energy intake and consumption. Parameters α and b are positive coefficients associated to energy intake and energy consumption for maintenance, respectively. Note that, while energy consumption is proportional to the whole volume of the organism (bV), energy intake depends only on a fraction of this volume (αV β ), provided by the metabolic scaling exponent β. This β, defined hereafter as the growth exponent, influences the shape of the tumor growth curve (Fig. 6).
While exponents β that characterize allometric growth of organisms are typically lower than 1 17,18 , it has been argued, and supported with data from different human cancers, that malignant tumors with heterogeneous clonal composition have exponents β > 1, what lead to a super-exponential growth. The reason is that in heterogeneous tumors there would be a range of phenotype/genotypes to choose from, leading to selection for more aggressive phenotypes/genotypes, and an increase of the growth exponents, which would manifest in the form of fasterthan-exponential explosive unbounded growth 8 . In this paper we will focus on scenarios of growing tumors for which the first term in Eq. (1) dominates over the second. In this context, it will be assumed that b ≃ 0, meaning that most of the metabolic requirements are routed towards biosynthesis rather than basal energy consumption. Moreover, our dataset having time intervals with three consecutive longitudinal measurements without treament (V 0 , V 1 , V 2 ) allows us to identify at most three parameters for each tumor (V 0 , α, β) but not more. When b is assumed to be zero, for a clonally homogeneous tumor with only a fraction of the cells proliferating due to necrosis, nutrient limitations, etc., exponents 0 < β < 1, correspond to growth slower than exponential, but with volumes still increasing with time.
Three patient groups were studied in a first batch of analyses. They included the cases of (i) growing untreated and (ii) relapsing post-radiotherapy BMs, and also (iii) growing BMs from patients under chemotherapy (CT) but with no specific treatment for the BMs. Patients in the last group included only drugs crossing the blood-brain-barrier (BBB) as described in 'Methods'. Figure 1a, b shows examples of BM longitudinal growth dynamics as observed in MRI studies. The growth exponent β governing the dynamics for each patient was obtained as described in the Methods section. The median value of the fitted individual exponents for untreated BMs (n = 10) was β = 1.59. This suggests substantial growth acceleration with β ≃ 3/2. It is interesting to note that this number differs from the value 5/4 obtained from metabolic scaling data of primary tumors 8 .
Median growth exponent for BMs growing under CT was β = 0.64 (n = 16). For BMs growing after completing radiation therapy (RT) we obtained β = 0.63 (n = 30). Finally, for those having received RT and under CT we got β = 0.65 (n = 40), where RT can stand for WBRT, SRS or a combination of both. Box plots are shown in Fig. 1c. Similar growth exponents were obtained after receiving different RT modalities (see Fig. 1d). Thus, volumetric growth of treated BMs had lower β exponents on average than those obtained for untreated BMs.
Animal models recapitulating BM's natural history display superexponential dynamics To investigate the growth patterns of untreated BMs in faithful animal models, experiments in mouse models were performed as described in Methods. For the experiments, the human lung adenocarcinoma brain tropic model H2030-BrM3 was used since it is injected into the hearts of nude mice and leads to the formation of brain metastases from systemically disseminated cancer cells. Thus, the cell line used recapitulates at least in part the complex sequences of transformations required for cells to metastatize.
Since more than three points were available in the mice dataset, the growth exponents were computed using a different fitting technique as described in 'Methods'. The median value of the individual exponents β was 1.44 (n = 20) (Fig. 1e, f).

Sensitivity analysis of the exponents calculation
Since three points per patient were available to obtain three parameters, the fitting could be very sensitive to small variations in the data. Those variations in volumetric data could be given either by the time between MRIs or by image segmentation, regardless of being performed by the same image expert, and revised by another expert and a radiologist. In order to asses the effect of small changes in measured volumes on the results, an analysis of sensitivity was performed as explained in 'Methods'.
The growth exponent values were consistent in 82% of the BMs used in the study when a random error was added. Once sensitive cases were excluded, the results were in line with those of the entire cohort of BMs (Fig. 2), ensuring the strength of the study.
Growth exponents best fitting the dynamics display superexponential growth for untreated patients and subexponential growth for treated ones An alternative approach to obtaining β was designed by looking for the value β * that provided the best fit for all patients in each subgroup: 1) untreated, 2) growing while in CT treatment (CT), and 3) relapsing post-RT -with or without CT− (RT*). The absolute errors weighted by volume, relative errors, were computed for each β * value and each group of patients.
The exponent best fitting the whole dataset of untreated BMs was β * = 1.5 (n = 10). For treated relapsing tumors, the best fit was obtained with β * = 0.51 (CT, n = 16) and β * = 0.71 (RT*, n = 70), showing again a slowing down of the post-treatment growth dynamics (Fig. 3b).   It's interesting to note that at β * = 1.5, a minimum appears in each dataset. Such minima are absolute for untreated lesions and relative in the case of treated ones, which indicates that after treatment, when the best growth exponent defining the group is lower than one, there is still a relevant tumoral component in some cases.

Evolutionary dynamics of tumor complexity
It is known that cancer treatments lead to a reduction of the clonal complexity during treatment, due to the selective pressures exerted by the therapies [19][20][21][22] . The successive overtaking of progressively fitter clonal populations is a key ingredient leading to the increase of the growth exponent β 8 . Hence, one would expect this reduction in heterogeneity to be reflected in the growth exponents of BMs as observed in our previous analysis.
The influence of different clonal compositions on the growth exponent of BMs was assessed in silico using an adapted version of the mesoscopic model developed in ref. 23 (as described in Methods). We explored computationally a simplified scenario were BMs were made of two predominant clonal populations, one being more aggressive than the other. The initially most abundant population proliferated and migrated at fixed rates, while the initially less abundant population had an advantage in both processes, and hence was assumed to be more aggressive, due to either mutational changes or irreversible phenotype changes, providing evolutionary benefits. In that way we accounted for an initial tumor heterogeneity, and the abundance of each population would change over time (due to one being fitter and more aggressive than the other, but less abundant), with the turning point taking place several months post SRS.
First, the growth exponent β of an untreated virtual BM was evaluated. To do so, we simulated starting with a small fraction (10%) of an aggressive population coexisting with a larger population of less aggressive cells. After a few months, the tumor was substantially enriched in silico in the most aggressive population (94% versus 6%) and the growth exponent β was found to be β = 1.53 Fig. 4a.
In a second set of computational experiments, we studied the effect of three scenarios after SRS treatment: a total depletion of the most aggressive subpopulation, a total depletion of the less aggressive subpopulation, and a balanced reduction in cell number of both subpopulations. In the first and second cases, that yield an homogeneous tumor after SRS, growth exponents were respectively β = 0.63 and β = 0.68, thus far from exponential growth (Fig. 4b, d). In the third case, we observe an heterogeneous tumor after SRS (a more realistic scenario), since we assumed that the most aggressive population was more sensitive to treatment, thus restoring equilibrium between the two subpopulations. In this scenario, the growth exponent found was β = 1.01, still far from super-exponential (Fig. 4c).
To study the influence of different advantages in proliferation and migration of the most aggressive population on the growth exponent β, simulations were performed varying the value of the advantage coefficients v div and v mig (see 'Methods'), and keeping the initial proportion of aggressive cells equal to 10%. We observed that the largest growth exponents β were obtained for combinations of slight advantages in both processes (β = 1.758 for v div = 1.25, v mig = 1.05; Fig. 5b) and that, for most combinations of coefficients v div and v mig , the resulting exponent β is lower than 1.
Therapy is most harmful to the most proliferative cells 24,25 . Hence, it is to be expected that the relative abundance of aggressive cells will decrease after treatment. To explore the influence of treatment on the growth exponent β, we repeated the previous analysis but this time the initial proportion of aggressive cells was modified to represent scenarios were the treatment was either very effective (initial proportion of aggressive cells is 1%, Fig. 5c) or not very effective (initial proportion of aggressive cells is 20%, Fig. 5a). The resulting β maps (each spot corresponds to a simulation performed with a fixed pair of Fig. 4 Simulations of longitudinal growth of heterogeneous BMs with two initial populations (turquoise: less aggressive, and ocher: more aggressive). a Pre-treatment and b-d post-treatment cases. The more aggressive population carries an advantage of 80% in proliferation speed and 92.5% in migration speed, compared to the less aggressive population. In a the BM is composed of 10% of more aggressive cells, and 90% of less aggressive cells. After eight months, the more aggressive population has overcome its counterpart, becoming dominant. Then, three different situations that can happen after treatment are illustrated: b the less aggressive population is completely removed from the tumor; c both populations remain in a balanced state, and d the more aggressive population is completely removed from the tumor. The betas were computed by choosing a random time point from each third of the total simulated time and are shown on each subplot. coefficients v div and v mig ) are qualitatively similar, although it can be observed that the maximum exponent β achieved decreases as the initial proportion of aggressive cells increases. Again, we observed that most combinations of advantage coefficients yield a β lower than 1. This result suggests that, regardless of the resulting proportion of aggressive cells after treatment, it is rare to observe a growth exponent β larger than 1 in treated BMs.

DISCUSSION
Evolution is one of the main driving forces of life on Earth and is behind the observed diversity at every level of biological organization. Evolutionary processes are used by cancers to survive within their hosts and escape from the pressures exerted by treatments. It is a remarkable fact that the growth laws of untreated human malignant cancers and their animal model counterparts display a signature of the evolutionary processes taking place behind the scenes, in the form of an exponent β > 1 in Eq. (1) 8 .
This study mined a dataset including more than a thousand BMs to test such a surprising result over a time scale of months, i.e. the time interval spanning three MRI studies (6-9 months). BMs have a background of heterogeneity that could provide the necessary substrate for evolutionary competitive dynamics to happen, leading to super-exponential growth of the tumor mass. Phylogenetic analyses have revealed that BM-competent clones genetically diverge from their primary tumors at a relatively early stage in lung adenocarcinoma patients 26 . Genomic analyses of solid tumors and matched BMs revealed significant genetic heterogeneity between primary lesions and BMs 27 , and the degree of genetic heterogeneity of BMs varied significantly among individuals with NSCLC, breast, and colorectal cancer [28][29][30] . In addition to the genetic heterogeneity, BMs have significant epigenetic variability 31,32 and there can be other phenotypebased mechanisms playing a role 33 .
Our results manifested both a macroscopic reflection of the evolutionary dynamics in the form of a β > 1 exponent for pretreatment longitudinal dynamics, and also the loss of biological richness experienced by BMs after therapy, which lead to substantially reduced exponents β ≈ 2/3 post-SRS. It has been hypothesized, using mathematical models, that treatment strategies in which an oncological "first strike" reduces the size and heterogeneity of the population, then followed by "second strikes" could lead to cancer extinction in metastatic disease 34,35 . Our results show the effectiveness of the first-line SRS approach in providing an ecoevolutionary first-strike strategy for BMs. Independently of the volumetric reduction observed, which ends up being marginally irrelevant if the tumor recurs, the substantial reduction of the growth exponent implies a direct effect on the tumor ecological complexity. In the case of BMs, "second-strike" strategies could be provided by targeted therapies with better penetration than classical drugs across the blood-brain barrier, many of which are under investigation 36,37 . It is also very interesting that information obtained from global macroscopic images could provide information on the underlying biological richness of these metastatic lesions. Thus β could be understood as an evolutionary exponent providing some information on the tumor heterogeneity.
It could be thought that the information about growth dynamics provided by the growth rate of a tumor (α in Eq. (1)) could be more important than that of the growth exponent β. As a toy example, let us consider two different BMs, one of them having a faster growth rate. For some values of the growth exponent β of the BM with the slower growth rate, it could happen that the "apparently slower" BM would end up surpassing (namely, growing faster in the long term) the "apparently faster" BM. Hence, the growth exponent β provides valuable information about the growth dynamics and the aggressiveness of BMs, independent of that provided by the growth rate. Moreover, this β also has practical implications in the management of BMs, as it can help in distinguishing between tumor recurrence and radiation necrosis. Radiation necrosis (RN) is a common side effect of irradiation that appears about a year after treatment and it can resemble true progression in appearance and clinical symptoms. It has been shown that growth exponents in the case of RN are larger than those observed not only in recurrent but also for untreated BMs 38 .
Regarding its prognostic value, one could think that the information provided by the β exponent may be biased by tumor size. An older BM may have had plenty of time to grow, so that the correlation between advanced time and BM size could be confounding the analysis of the growth exponent β. However, it should be pointed out that this exponent is an intrinsic property of the growth curve, so that it remains the same during the whole lifespan of the BM, and independently of the times of observation (three at least, in order to be able to fit Eq. (3) as described in Methods). As an example, such correlation between advanced time and larger size would still be observed for BMs with β = 1, which would be showing no acceleration nor explosive growth. Hence, β prognostic information is independent of tumor size.
It may seem naively counterintuitive that BMs subject to different treatment modalities (SRS, WBRT, CT) led to similar growth exponents β, since SRS is known to be substantially more effective than WBRT or CT. Our retrospective study focused only on BMs growing after (RT) or under (CT) treatment, but of course there would be many BMs with complete response, e.g. to SRS, Fig. 5 Growth exponents β obtained from a parameter sweep varying the advantages in migration and proliferation. Simulations were carried out for different combinations of coefficients v div (advantage in proliferation; values explored range from 1.11 to 5) and v mig (advantage in migration; values explored range from 1.05 to 5). a Poor-effectiveness treatment case (proportion of aggressive cells after treatment is 20%). The largest exponent β obtained is equal to 1.615, for a v div = 1.33 and a v mig = 1.05. b Medium-effectiveness treatment case (proportion of aggressive cells after treatment is 10%). The largest exponent β achieved was 1.758, for a v div = 1.25 and a v mig = 1.05. c High-effectiveness treatment case (proportion of aggressive cells after treatment is 1%). The largest exponent β achieved was equal to 1.99, for a v div = 1.1765 and a v mig = 1.05. Gray lines correspond to β = 1.
that were not included here. It is also relevant to emphasize that the growth exponent cannot be directly interpreted as a growth rate, e.g. the speed of volumetric growth, but as a measure of the shape of the tumor growth curve.
An intriguing result of our study was the fact that the growth exponent for untreated BMs when fitted together was close to 3/2. It has recently been found for different primary tumors, lung (both adenocarcinoma and squamous cell), breast, colorectal, glioma, and head and neck cancer, that metabolic scaling exponents are close to 5/4 8 . Following the classical reasoning of West et al. 39 , one would expect metabolic exponents to be the same as growth exponents. This raises the interesting question of what the metabolic scaling of BMs would be on diagnosis, since the datasets of 8 did not include that condition. Would it be 3/2, raising the question of why BMs have a different metabolic scaling than other cancers? Or would it be 5/4, raising the question of why there should be a mismatch between metabolic and growth exponents in BMs? Data on metabolic scaling of BMs would be necessary to answer that question.
The in silico observation that the largest exponents β were achieved for combinations of slight advantages in both processes, division and migration, points out that a great advantage will not lead to a β > 1, as the overtaking time will be reached more quickly by the aggressive population. Another observation is that β decreases faster with changes along the v mig axis, suggesting that advantages in migration may bring more competitive advantages in division.
Several authors have wondered whether the interval between SRS planning and treatment is accurate [40][41][42][43] . In order to study this, they compared volumes from diagnostic imaging and radiosurgery planning MRI and extrapolated the growth linearly to the day of SRS, since only two time points were available. Progression between diagnosis and SRS is common, and some suggest that a mathematical model would be useful to individualize treatments. Our model based on three time points shows that the growth velocity of untreated BMs increases with time, thus pointing out an inadequate prediction of the tumor volume on treatment day and a substantial benefit of reducing the interval between SRS planning and treatment.
The main strength of our study lies in the utilization of 96 brain metastases (BMs) from a substantial dataset comprising 1133 BMs treated with stereotactic radiosurgery (SRS), characterized by highresolution data and adherence to the guidelines for BM clinical studies 44 . The selected subset of BMs was suitable for computing β as they exhibited sustained growth across three sequential imaging studies, without undergoing additional therapeutic interventions during that time window. Additionally, each lesion was meticulously segmented by the same expert and verified by a radiologist. Moreover, our study adopted a multicenter approach, encompassing BMs from five different hospitals.
However, a key limitation of our study is that due to the size of the final dataset and the predominant representation of lung cancer cases, we were unable to perform separate analyses by primary histologies. It would be of great interest to investigate whether the conclusions drawn from our study are dependent on the specific type of primary cancer considered.
In summary, we studied a large BM dataset and unveiled a continuous acceleration of growth in the case of untreated lesions, due to evolutionary dynamics sustained between different tumor subpopulations, as validated by in silico simulations using a stochastic discrete mesoscopic model. Results for mice data were in line with that. Recurrent BMs after treatment displayed slower growth, compatible with treatment-mediated reduction of tumor heterogeneity. As a result, we have highlighted the predictive value of a macroscopic variable, the growth exponent, which can be used to obtain information on the microscopic status of the tumor.

METHODS AND MATERIALS Patients
Patients included were all participants in the study MetMath (Metastasis and Mathematics), a retrospective, multicenter, nonrandomized study approved by five hospitals. All patients were diagnosed with BM in the period 2007-2021 and followed up with MRI according to standard clinical practice. A total of 354 patients who received SRS at any time during the evolution of the disease, and with full longitudinal follow-up, were reviewed in the study, including 1133 BMs. Primary tumor histologies were mainly nonsmall-cell lung cancer (NSCLC), breast cancer, melanoma and SCLC.
For the study of longitudinal volumetric dynamics, BMs were selected from the MetMath dataset on the basis of several inclusion criteria. First of all a minimum of three consecutive imaging studies, including a volumetric contrast-enhanced (CE) T1-weighted MRI sequence (slice thickness ≤2.00 mm, no gap) with no substantial imaging artifacts, at different time points, were required in order to allow for reliable lesion volume calculations. Secondly, an increase in tumor volume at each of the three time points was required, since it was desired to study the growth of either untreated or recurrent tumors. Next, only time points without previous SRS/WBRT treatments (for untreated cases) or with SRS/WBRT treatments received more than four months before the first imaging study were considered (treated cases), in order to exclude the potential confounding effect of acute inflammatory responses seen in some patients in the first MRI were also excluded. Not all chemotherapies (CTs) are able to cross the blood-brain-barrier (BBB), including pertuzumab or trastuzumab 45,46 . Such CTs target primary tumors but are known to have no effect on metastatic lesions due to the protecting effect of the BBB. Lesions that are not irradiated and with drugs known as not able to cross the BBB, were consider as untreated. Each BM in our dataset was carefully revised to determine whether or not it satisfied the inclusion criteria. To do so lesion segmentation was necessary in many cases to assess its growth dynamics. Finally, 96 BMs from 69 patients were included in the study. Of these, 10 were untreated, 16 had received chemotherapy (CT), 30 WBRT or SRS (radiation therapy) and 40 received both treatments. A summary of patient characteristics used for the study is provided in Table 1.

Imaging and follow-up
The volumetric contrast-enhanced T1-weighted MR imaging sequence used to delineate the BMs and compute their volumes was gradient echo using 3D spoiled gradient-recalled echo or 3D fast-field echo after intravenous administration of a single dose of gadobenate dimeglumine (0.10 mmol/kg) with a 6-to 8-minute delay. MRI studies were performed in the axial or sagittal plane with a 1.0 T (n = 5), 1.5 T (n = 357) or 3.0 T (n = 55) MR imaging unit. Imaging parameters were no gap, slice thickness of 0.52-2.0 mm (mean 1.3 mm), 0.4-1.1 mm (mean 0.5 mm) pixelspacing and 0.4-2.0 mm sacing between slices (mean 1.0 mm).
Typical time spacing between MRI studies for BM follow-up was about 3 months for the institutions participating in the study. In our dataset, the median time between the first two MRI studies was 3.04 months while for the second it was 2.64 months. . Each BM lesion was automatically delineated using a gray-level threshold chosen to identify the CE tumor volume. Segmentations were then corrected manually, slice by slice, using an in-house software as described in ref. 47 . Necrotic tissue was defined as hypointense tumor regions inside CE tissue. CE and necrotic areas of the lesions were reconstructed, the tumor interfaces rendered in 3D. Tumor volume was computed as the volume within the surface delimiting CE areas.

Experiments in animal models
The human lung adenocarcinoma brain tropic model H2030-BrM3 (abbreviated as H2030-BrM) 48 was injected into the hearts of nude mice to induce the formation of brain metastasis from systemically disseminated cancer cells. H2030-BrM was cultured in an RPMI1640 medium supplemented with 10% FBS, 2 mM l-glutamine, 100 IU ml −1 penicillin-streptomycin and 1 mg ml −1 amphotericin B. Brain colonization and growth of metastasis were traced using noninvasive bioluminescence imaging, as BrM cells express luciferase. On administration of the substrate D-luciferin, bioluminescence generated by cancer cells was measured over the course of the disease. The increase in photon flux values is a well-established correlate of tumor growth in vivo 48 . The experiments were performed in accordance with a protocol approved by the Centro Nacional de Investigaciones Oncológicas (CNIO), the Instituto de Salud Carlos III and the Comunidad de Madrid Institutional Animal Care and Use Committee. Athymic nu/nu mice (Harlan) aged 4-6 weeks were used. Brain colonization assays were performed by the injection into the left ventricle of 100 μl of PBS containing 100,000 cancer cells. Mice anaesthetized with isofluorane were injected retro-orbitally with D-luciferin (150 mg kg −1 ) and imaged with an IVIS Xenogen machine (Caliper Life Sciences). A bioluminescence analysis was performed using Living Image software (v.3).

Ethical approval
We have complied with all relevant ethical regulations. Human data were obtained in the framework of the study MetMath (Metastasis and Mathematics), a retrospective, multicenter, nonrandomized study approved by the corresponding institutional review boards: Fundación Instituto Valenciano de Oncología, Hospital Universitario HM Sanchinarro, Hospital Regional Universitario de Málaga, MD Anderson Cancer Center and Hospital Universitario de Salamanca.
Animal care and experimental procedures were performed in accordance to the European Union and National guidelines for the use of animals in research and in accordance with a protocol approved by the CNIO and Comunidad de Madrid Institutional Animal Care and Use Committee (H2030-BrM3 cells).

Von Bertalanffy growth model
Solving (1), with b = 0 leads to where V(t) is the volume (as a function of time t) of a BM, V 0 is the volume obtained from the first segmentation of a BM, α is the energy intake coefficient, and β is the metabolic scaling exponent (as in Eq. (1)). Since there is information about the dynamics at three time points (t 0 , V 0 ), (t 1 , V 1 ) and (t 2 , V 2 ) obtained by image segmentation, the two parameters α and β can be completely determined by evaluating (2) at the times t 1 , t 2 , giving Equation (3) is an algebraic equation for β that was solved using the MATLAB function fzero (which returns the root of a nonlinear function) for each set of known values V 0 , V 1 , V 2 , t 0 , t 1 , t 2 . A sensitivity analysis was performed to ensure the robustness of the method for the computation of β.
The growth exponent β provides information on the shape of the tumor growth curve, which cannot be directly interpreted as a growth rate, e.g. the speed of volumetric growth. Figure 6a shows examples fitting the same pair of volumes and time points (initial and final), with different values of β while subfigure 6b shows that for a fixed value of α, beta represent the aggressiveness of the tumor. Figure 6c, d, illustrates that any fixed value of β (chosen there arbitrarily as β = 1, i.e. exponential growth and β = 5), is compatible with different speeds of volumetric growth.

Sensitivity analysis of the exponents calculation
To ensure the robustness of the growth exponent β computation, some well-known growth curves were studied: exponential, cubic, Gompertz and logistic. For each type of growth, 100 sets of points B. Ocaña-Tienda et al.
(t 1 , t 2 , t 3 ), (V 1 , V 2 , V 3 ) were randomly chosen and β was computed. For the exponential and cubic growths, the same value was obtained independently to the chosen points ( Fig. 7b1-b2). In the case of Gompertz and logistic growths a set of β values was obtained (Fig. 7b3-b4). When adding a random error from -20% to +20% to the curves, a wider range of values for β appears, but such values are located around the values obtained without error (Fig. 7c).
Longitudinal growth analysis: individual growth exponents for mice In animals, more than three volumetric points were available and a different method was used to calculate the growth exponent for each mouse. Measurements at 7, 14, 18, 21, 25 and 28 days were usable from 24 mice, however at 25 and 28 days, the measured volume was close to the total brain volume and it was assumed that growth would be affected by limitations of space. Thus, to avoid confounding effects coming from mechanical constraints, the latter two time points were excluded from the analysis. A discretization of Eq. (1) was performed, and taking logarithms on both sides the slope of the straight line that best fits all the points for each mouse corresponds to the growth exponent β. Four mice were excluded because of volume decrease, leading to the use of 20 mice for β computation.

Sensitivity to small changes in volume when computing beta
In order to test the sensitivity to little fluctuations in the data, a random error smaller or equal to ±5% was added to each volume. For each BM, the procedure was carried out 200 times in order to compute β * for each set of random errors. It was imposed that the average of the 200 calculated β * have a difference less than 0.5 from the computed β for the measure volumes, that is to say, jβ Ã av À βj < 0:5.

Longitudinal growth exponents: Group calculations
An iterative method was used to automatically compute the optimum β * , that is, the one giving the lowest relative error to the segmented volumes for each group. A sweep was performed on β 0 ¼ ½0; 3 with 300 steps and on V 0 0 ¼ ½0; V 1 and V 0 1 ¼ ½V 0 ; V 2 , with 500 steps for each. To each β 0 and each pair ðV 0 0 ; V 0 1 Þ there corresponds a single value of V 0 2 from Eq. (3). Then, the sum of relative errors of the three segmented volumes for every BM (V 0 , V 1 , V 2 ) was computed using the formula: the smallest is retained, that is to say, the combination of ðV 0 0 ; V 0 1 ; V 0 2 Þ that best fits (V 0 , V 1 , V 2 ), given each β 0 (pseudocolor plots in Fig. 3a). Finally, the β * value for which the sum of all errors (Fig. 3b) is minimum and therefore corresponds to the best fit for the whole subgroup, is observed.
Stochastic discrete mesoscopic simulator of BMs longitudinal dynamics and response to treatments To illustrate how treatment could modulate heterogeneity and influence the β values obtained from BM longitudinal growth data, Fig. 6 Interpretation of the growth exponent β. a Growth behavior for several values of the growth exponent β when fixing initial and final volumes and times. The growth exponent gives information about the shape of the curves. b Growth behavior for the same values of the growth exponent β than in a when fixing initial volumes and α. c Growth behavior when fixing β = 1 and the initial volume, showing different growth curves with different α but the same exponent. d Growth behavior when fixing β = 5 and the initial volume, showing different growth curves with different values of α.
as well as to provide an in silico testbed capable of simulating the growth and evolution of BMs, an adapted version of the mesoscopic model developed in ref. 23 was used. The mesoscopic model is a discrete stochastic simulator of tumor growth, that features clonal populations as the basic agent, instead of individual cells. The spatial domain is divided into voxels, a 3D generalization of a bi-dimensional compartment (commonly used in medical imaging). Each voxel has a given carrying capacity K, so it can harbor as much as K cells, that may belong to different clonal populations. It is assumed that cells belonging to the same clonal population and exposed to the same stimuli (microenvironment, surrounding cell density, nutrient/oxygen availability) will behave in the same way. Therefore, instead of assessing cell-bycell the outcome of any dichotomous success/failure process that an individual cell may perform (such as division or death), the number of cells in each clonal population that succeed in performing such processes is assessed in a single step.
If the outcome of each individual process i performed by a cell is considered as a random variable following a Bernoulli distribution X i~B ernoulli(P i ), with a probability P i associated to that process, then the joint outcome of an entire clonal population of N identical cells performing such process can be considered as a random variable following a binomial distribution X 0 i $ B ðN; P i Þ. In this way, by knowing the number of cells N attempting to perform a process i, and the probability P i of successfully performing that process, updating the number of cells in a clonal population is done by sampling the binomial distribution associated with that process.
The processes considered in this version of the model are division, death and migration. Note that mutations or phenotypic transitions are disregarded, as cells are not allowed to change from a given population to another. This assumption is grounded in the short time scale of evolution of BMs. At each time step, the cell number in each voxel is updated in a synchronous way (using swapping matrices) by random sampling the binomial distributions corresponding to each process, voxel and clonal population. The probabilities associated with these binomial distributions were defined in the same way as in the original work 23 . The version of the model used in this paper is adapted to simulate BMs. It considers only two clonal populations, with different traits and characteristic rates. The rates of division, death and migration were fitted (using ABC rejection algorithm 49 ) to mimic the lifespan and volume dynamics of an actual BM.
One of these clonal populations is more aggressive than the other; this is implemented by faster division and migration rates. Namely, the advantage associated with the division process will be v div ∈ [1, ∞), while the advantage associated with the migration process will be v mig ∈ [1, ∞). These coefficients represent the ratio between the characteristic rate (inverse of time) at which the aggressive cells carry out the considered process, versus the characteristic rate of the other cell type. Hence, the greater their value, the greater the advantage associated with the aggressive population. As an example, a v div = 2 would mean that the rate of division of the most aggressive population is twice the rate of division of the less aggressive population; in other words, the division time of the most aggressive population is half the division time of the less aggressive population. The range of parameter values used in the model to perform simulations can be seen in Table 2. A total of 720 simulations were performed for the parameter exploration performed in Fig. 5.
The adapted version of the mesoscopic model for BMs was coded in Julia (v. 1.1.1). Data processing and visualization of simulation files was performed in Matlab (v. 2018b). The original code for the mesoscopic model has been made available in the repository https://github.com/JuanJS117/MesoscopicModel. Fig. 7 Growth exponent β for different types of growth: exponential, cubic, Gompertz and logistic. a Growth function in blue with a ±20 % error in green. b Computed β from Eq (1) for the different growth without error. c Computed β when a random error is taken into account.

Statistical analyses
Statistical analyses were performed using the MATLAB software, and also SPSS (Statistical package for the Social Sciences, v24.00 IBM) software. The normality of the variables was assessed via the Kolmogorov-Smirnov test. The Kruskal-Wallis test was conducted with adjustment for multiple comparisons, to determine statistically significant differences for non-parametric data (the scaling law growth factor, β). P-values smaller than 0.05 were considered to be statistically significant.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.

DATA AVAILABILITY
All data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.